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A comprehensive microscopic dynamical theory is presented for the description of quantum fluids as they 
transform into glasses. The theory is based on a quantum extension of mode-coupling theory. Novel effects 
are predicted, such as reentrant behavior of dynamical relaxation times. These predictions are supported by 
path integral ring polymer molecular dynamics simulations. The simulations provide detailed insight into 
the factors that govern slow dynamics in glassy quantum fluids. Connection to other recent work on both 
quantum glasses as well as quantum optimization problems is presented. 



I. INTRODUCTION 

Understanding the fundamental causes of the dramatic 
slowdown of dynamics when a liquid transforms into 
a glass is still a subject of great debateP^ Essentially 
all discussion of the glass transition has focused on the 
strictly classical regime of liquid state behavior, namely 
where the de Broglie wave length is significantly smaller 
than the particle size. Given that nearly all known 
glass forming liquids fall well within this regimeP it is 
clear that the classical approximation is generally jus- 
tified. However there are several interesting and impor- 
tant examples where quantum fluctuations and glassiness 
coexistP-^ In such cases, which range from the behavior 
of superfluid helium under high pressure to the phase di- 
agram of quantum random optimization problems, the 
interplay between quantum mechanics and the otherwise 
classical fluctuations that lead to vitrification can be ex- 
pected to produce qualitatively novel physical behaviorP 

The theoretical investigation of quantum glasses has 
increased in recent years. Studies ranging from the in- 
vestigation of quantum effects in so-called stripe glasses, 9 
quantum spin-glasses^HUl anc l lattice models that mimic 
the properties of superfluid and supersolid heliurrf^ have 
been presented. In this work we instead focus on "realis- 
tic" off-lattice quantum fluids. The microscopic detail of 
our study necessitates the use of approximations, such as 
mode-coupling theory (MCT)P^ and ring-polymer molec- 
ular dynamics (RPMD) P3 that are less well-justified then 
the methods employed in the studies of the model sys- 
tems mentioned above. On the other hand, the ap- 
proaches used here have lead to a host of non-trivial pre- 
dictions both for classical glass-forming liquidglS] as well 
as a variety of quantum liquid-state phenomena.^ We 
thus expect that the predictions made in this work to be 
at least of qualitative accuracy. 

The work presented here builds on our earlier report 
of several novel effects that arise when glassy dynam- 
ics occurs in the quantum regimeP In particular, both 



RPMD and the quantum version of mode-coupling the- 
ory (QMCT) indicate that the dynamical phase diagram 
of glassy quantum fluids is reentrant. As a consequence, 
hard-sphere quantum liquids may be forced deeper into 
the glass "phase" at fixed volume fraction as quantum 
fluctuations increase. This counterintuitive finding has 
implications not only for liquid-state systems such as su- 
perfluid helium under pressure, but for a broad class of 
quantum optimization problems as well. 

In comparison to our earlier paper P the work presented 
here provides complete details for both the QMCT and 
the quantum integral equations needed for generating the 
required structural input. In addition, we give a far more 
extensive interpretation of the results, largely afforded 
by our RPMD simulations. Lastly, we discuss in greater 
detail the connection of our results to related theoretical 
work. 

The paper is organized as follows: In Sec. [IT] we provide 
the details of the QMCT, including a description of the 
equations for the density correlator and the mode cou- 
pling approximations. In addition, we discuss the high 
and low temperature limit of the QMCT and derive equa- 
tions for the nonergodic parameter used to determine the 
liquid-glass line. In Sec. |III| we describe the quantum in- 
tegral equation theory used to obtain the the static in- 
put required by QMCT. Sec.|TV]is devoted to the RPMD 
method. Results and discussions are presented in Sec.|VI| 



Finally, in Sec. VII we conclude. 



II. A SELF-CONSISTENT QUANTUM MODE-COUPLING 
THEORY 



The general quantity of interest is the Kubo trans- 
forrrPH of the time correlation of the collective density 
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operator, p q = J2a=i e jqfc> , given by 
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with a time evolution described by the exac^ quantum 
generalized Langevin equation (QGLEpSl 

4> q (t)+tf q <f> g (t)+ [ drM q (T)4> q (t-r) = 0, (2) 
Jo 

In the above, we have used the notion that r a stands for 
the position vector operator of particle a with a conjugate 
momentum p Q and mass m, N is the total number of 
particles, /? = is the inverse temperature and (• ■ • ) in 
Eq. ([TJ) denotes a quantum mechanical ensemble average. 
The frequency and memory terms are given by: 



n 2 = ^ 

q mpcj> q (Q) 



and 



M q {t) = 
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respectively, with C = jjr [if, • • • ] being the Liouvillian 
and £ = QiQ\HQ\Qi- To derive the above equations we 
have defined two projec tion operators (first and second 
order, respectively) 22 * 23 * 



and 



Pi= \p q )^\0)( Pq \ (5) 



P2 = \QiC Pq ) {QiCp q \Q x Cp q ) 1 (Qi£p q \ (6) 



with Qi = 1 — Pi and Q2 = I — P2. <f> q (0) is the zero time 
value of (j> q (t) and can be approximated bjBH -^JLa 

where 5 g is the static structure factor, An(w) = n(ui) — 
-311^ — r is the Bose distribution fulle- 



rs— a>) and n{uS) 
tion at temperature T. 



A. Quantum Mode-Coupling Approach 

We employ a quantum mode-coupling approach re- 
cently described by us for quantum liquids^ to obtain 
the memory kernel described by Eq. Q. This approach 
is based on our early work to describe density fluctu- 
ations and transport in quantum liquids such as liq- 
uid para- hydrogen , ori/io-deuterium, and normal liquid 
heliumPSES The basic idea behind this approach is that 
the random force projected correlation function, which 
determines the memory kernel for the intermediate scat- 
tering function, decays at intermediate and long times 
predominantly into modes which are associated with 



quasi-conserved dynamical variables. It is reasonable to 
assume that the decay of the memory kernel at long times 
will be governed by those modes that have the longest re- 
laxation time. Thus, the first approximation made by the 
QMCT is to replace the projected time evolution opera- 
tor, e iCt , by its p roje ction onto the subspace spanned by 
these slow modesPS The second approximation involves 
the factorization of four-point density correlations into a 
product of two-point density correlation.^ 

Following the derivation outlined by Gotze and Liickc 
(GL) for zero temperature pUM the memory kernel 
at finite temperature (in frequency space), M q (uj) — 
dte~ lu)t M{q,t)), can be approximated by 
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where n is the number density, 4> q (oj) = J^*L dte lu)t (f) q (t) 
is the Fourier transform of the Kubo transform of the 
intermediate scattering function and 

T(wi, uj 2 ) = n(— a>i)n(— W2) - n(wi)n(w 2 ). (8) 

The vertex, V q ^, q ~k, is formally given by 

N q -k,kV q ,k,q~k = (Q£ 2 pq\PkPq-k) (9) 
= (£ 2 Pq\PkPq-k) - Q q (pq\pkPq-k) , 

with the normalization approximated by 

Nq-k,k = (p q Pq-k\pq-kPk) (0) (10) 
[°° du [°° djj 1 

Hp I — I — — T(u),l>-w) 



„ 7T J-oo 4w 

xuj'(lo - uj')4> q - k (uj')4> k (uj - u'), 

consistent with the spirit of QMCT where four-point den- 
sity correlations are factorized into a product of two-point 
density correlations P3 

B. The vertex 



The vertex in Eq. ( 10 1 is difficult to compute since it 



involves three-point Kubo density correlations. A com- 
mon approach taken by classical mode-coupling theory 
(CMCT) is based on a convolution approximation. 35 For 
the Kubo transform quantum case, a convolution-like 
approach is not unique. The approach we adopt here 
is based on an ex tension of the work of GL to finite 
temperatures) 2 ^ 3 ^ In this work, a dynamical approxima- 
tion is made to remove the dependence on Kubo trans- 
formed structure factor in the vertex. The assumption 
behind this approximation is that the major contribution 
to the vertex and its normalization comes from a charac- 
teristic frequency of the system. Thus, we approximate 
4> q (uj) within the vertex by 
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which satisfies the known sum rule dui<f) q (oj) = <fi q (0). 
Inserting this approximation for <p q (oj) into the expres- 
sion for N q _k,k given by Eq. (11) yields: 
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2S q — k S k 



hf3An({l q _ k )An(Vl k ) 



K(n q _ k ,n k ), (12) 



where 

v(Vl x r(O g -fc,^ fc ) T(-fl q - k ,tt k ) 

"g-fc + "fe "g-fc ~ "As 

For V q ^ k _ q - k we use the exact relations^ 

(L 2 p q \p k p q ^k) = — (g • fc^-fe + q- (q- k)Sk) ■ (14) 
v ' mp 

and the convolution approximation 

(Pq(t),PkPq-k) ~ S q S k S q ^ k (15) 

to obtain the approximation to the vertex: 
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where 
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f2 g Art(f2 fc + fi g _ fe ) - (fifc + f2g_fe)An(Q g ) 



(17) 



The above expressions close the equation of motion 
(Eq. Q) and require only the static structure factor to 
produce a full approximation to the time dependence of 
the quantum density-density time autocorrelation func- 
tion. 



duces to: 
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- uj')4> q ^ k {uj')4> k (uj - u)'), 



with 



C. High and Low temperature limits 

It may be shown that the above equations reduce to the 
venerable classical mode-coupling equations in the high 
temperature limit and to the GL theory as T — > 0. The 
latter theory produces a representation of the dispersion 
of superfluid helium that is at least as accurate as the 
Feynman-Cohen (FC) theory^ at low values of q and 
exhibits Pitaevskii-bending of the spectrum at high q, 
unlike the FC theory. In particular at high T, 
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d 3 k (q ■ kc k 
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+ q- (q- k)c q - k ) g _ k (t)0 fc (i), 



where c q = ^ ( 1 



is the direct correlation function. 



In addition, 4> q (t) reduces to the classical intermediate 
scattering function, F(q, t) as j3 — > 0. This is recognized 
as the CMCT memory function!^ 

At T — > the equation for the memory function re- 



lim V q k n- 



hn 
2m 

x (q ■ kc k 



(0J k + Wq-fe + ^q) 

• (<? - k)c q -k) 



(20) 



which are the T — > equations for quantum densi ty flu c- 
tuations in superfluid helium first derived by GL) 22 * 34 ! In 
the above, ui q — . We note in passing that the term 
(3 2 appearing in Eq. ( 20 I (and not in the derivation of GL) 



arises from our definition of the Kubo transform (Eq. [TJ , 
which includes a j^, while that of GL does not. Care 
must be taken applying the Kubo transform as T — > 0. 

In the T — > case, the entire structure of the memory 
function differs greatly from that of its high tempera- 
ture counterpart and the convolution structure is lost. 



Eqs.(20l and (21) do not imply a memory function that 
is a product of correlators at identical times. This is a 
consequence of the quantum fluctuation-dissipation the- 
orem (QFDT) that must be satisfied. At T -> the 
function T(iv q ,uj k ) becomes proportional to the difference 
of a product of step-functions in frequency, dramatically 
altering the structure of the theory. This distinction be- 
tween the low and high temperature limits has important 
consequences, as discussed below. 
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D. Nonergodic parameter 



is simple to show that f q must satisfy the equation: 



.38 



The nonergodic parameter, 
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HAn(Q, g )(l q 
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4> q (t-KX>), (21) 



is often used to describe the ergodic to nonergodic tran- 
sition as the liquid is cooled down to the mode-coupling 
critical temperature T c . Above T c one finds a single so- 
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The above equation for the nonergodic parameter re- 
flects the structure of the QGLE (Eq. Q), and thus, 
is valid both in the classical and quantal limits. In 
the former, the long time limit of the memory kernel 
is given by M q (t -> oo) « J d s k?l fc __*/,-*/* 



With V q,k,q-k 



Sq-kSk {q 



167r 3 mg 2 

kc k + q 



(q - k)c q - k ) ■ The 



quantum case is a bit more complicated since the struc- 
ture of the memory kernel is quite different and involves 



lution where f q = for all values of q, while at T c the a convolution of products of <p q (uj). The derivation for 
nonergodic parameter acquires a finite value f q > OP^It M q (t — > oo) is thus, based on the following expansion: 



-TV, uj - a/V (w - u') = — + ^uj'(uj - u') - (uj' 2 - u/(u - J) + (u- uj') 2 ) + 

uj V ' ' V J [3h 12 y ' 720 V 1 ' V M 
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(23) 



((w - uo'f - uj'(uj - uj') 3 + uj' 2 (uj - uj') 2 - uj' 3 {uj - <J) + uj' 4 ) + 0(/3 7 ). 



Inserting this into the memory kernel (Eq. ([8|) and keep- 
ing the first two terms only, we obtain: 

~ , , hm[3 2 f d 3 k , f°° . ( 1 

Y^w'(w - u>') H ^ <j) q - k {io')4> k (uj - uj'). 

In the time domain, this translates to: 



M ? (t) 



hm/3 2 



d 3 k 
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where the dot denotes a time derivative, i.e., 4>k (t) 



d ^ fc W The other terms in the expansion of Eq. (|23|) that 



<it. 



have been omitted give rise to terms of the form 

E^M-feW^W, (26) 

where a,(/3) are related to the expansion coefficients of 

iT(a/,w - w')fa/(w - a/) and ^(i) = i s the j's 

time derivative of 4> k (t). 

The long time limit of the Eq. (|26[) is now given by: 



M q (t -t OO) 



rn/3 f d 3 k 



V 2 

2q 2 n J (2tt) 3 r «' fc ^- fe 
X(j) q ^ k (t -> oo)cj>k{t -> oo), 



(27) 



where all the time derivatives vanish as i — >■ oo even when 
0fe(f — > oo) decays to a constant. Finally, we can rewrite 



the above in terms of the nonergodic parameter: 



M q (t -> oo) 



to/3 f d 3 k 



2q 2 n J (2tt) 3 

^V q 2 k ^ k ^ k (0)M0)f q -kfk- 



(28) 



The above expression is strictly valid at T — > but not 
at T = 0, since the expansion given by Eq. (23) is not 
valid at T — 0. The final result is similar to the classical 
equation, however the vertex is given by the full quantum 
mechanical expression of Eq. ( 16 I. 



III. QUANTUM INTEGRAL EQUATION THEORY 



The QMCT requires as input the static structure fac- 
tor, S q and its Kubo transform </>«{0). Here, instead of 
using PIMC to generate this input p2 we refer to a quan- 
tum integral equation approach, that i s bas ed on the 
early work of Chandler and Richardson! 40 ! 41 ! We begin 
with the Ornstein-Zernike relation applicable to quan- 
tum liquids. The quantum system composed of N parti- 
cles can be mapped on a classical system consisting of N 
ring polymers, each polymer being composed of P beads. 
Then, we can write the matrix RISM (reference interac- 
tion site modeP^^) equation for the classical isomorphic 
system by: 

h(\r — r'|) = w * c*w(\r — r'Q + nw * c* h(\r — r'|), (29) 

where * denotes a convolution integral and as before, 
n is the number density. In the above equation, h(r), 
w(r), and c(r) are the total correlation function, the self 
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correlation function, and direct correlation function, re- 
spectively, defined by: 



h{r) = 
w(r) — 
c(r) = 







hp 
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dXh(r,X) 
d\w(r, A) 
dXc(r, A), 



(30) 



and h(r, A), w(r, A), and c(r, A) are the imaginary time 
total, self, and direct correlation functions, respectively. 
In the classical limit Eq. ( 29 1 reduces to the classical 



Ornstein-Zernike equation with w(r) = 1. In what fol- 
lows, we will use the notation w q (X) for the Fourier trans- 
form of w(r, A), and similarly for c,(A) and h q (X): 

r hj3 

h„ 



1 

W Jo 
HP J 

i 

hp 



dXh q (X) 

dXWq(X) 

dXc q (X). 



(31) 



To proceed, we refer to the mean-pair interaction ap- 
proximations along with the quadratic reference action^ 
and rewrite: 



where 



w q (X) = exp{—qR(X)} 



i#r\\ _ >^ 1 ~ cos(0/A) 



(32) 
(33) 



to is the particle mass, Vtj — 2%j/hp is the Matsubara 
frequency and otj is given by: 



1 

6tt 2 hp 



h/3 



dq / dXq v q (l — cos(fljX)w q (X). 
Jo 

In the above the solvent induced self-interaction is given 
by: 



-Cq(nw q + n 2 h q ). 



In order to close the quantum Ornstein-Zernike equa- 
tions, which in g-space can be written as: 



WqCqWq + UWqCqhq, 



(36) 



we use the Percus-Yevick (PY) closure of the form (in 
r-space) : 



c(r) = (h(r) + c(r) + l)(exp(-/3u(r)) - 1), 



(37) 



where v(r) is the pair interaction between two particles. 
The static structure factor and its Kubo transform are 
then given by: 

nh q (38) 



S q = 1 



0,(0) = w q + nh q . 



In all the applications reported below we have used the 
approximate relation for 0,(0) ~ pnAn(n )n * 



IV. RING POLYMER MOLECULAR DYNAMICS 

The RPMD approach to quantum dynamics provides 
an approximation to quantum mechanical Kubo trans- 
formed correlation functions by using a classical evolu- 
tion of the imaginary time paths 19 . Consider a multidi- 
mensional system of N distinguishable particles with a 
Hamiltonian of the form, 



H 



N 

E 



2to q 



y(ri,...,rjv), 



(39) 



where, r a and p a are the positions and momenta of the 
particles and V(ri, . . . , r/v) is the potential energy of the 
system. The RPMD approximation to the canonical cor- 
relation function, CAE,{t), for position dependent opera- 
tors A(r) and B(v) is, 



CAB{t) 



1 



-73JVP, 



j3NP„ 



(2irh) 3NP Z P 
e- 0pHp ^Ap(r)B P (r t ) 



where 

Z P = 



{2nh) 



3NP 



j3NP r 



d 3NP r e -fc>*p(P.*) 



(40) 



(41) 



and Pp = PIP. Hp(p,r) is the classical Hamiltonian of 
the N particle P bead ring polymers with the external 
potential of V(ri, . . . , r/v) acting on each bead, 



N P /, (fe)s 2 i 
a=l k=l \ 



■ a / 



■ • i 1 jv ; 



(42) 



fc=i 



(35) where ujp = l/ph and the cyclic boundary condition 



ri P+1 ' 1 = r$ applies. The time-evolved coordinates 



r t = r t (p,r) in Eq. (40 1 are obtained from the classi- 
cal dynamics generated from this Hamiltonian and the 
operators Ap(r) and Bp(r t ) are evaluated by averaging 
over the beads of the ring polymer at times and t re- 
spectively, 



A P (v) 



k=l 



1 



k=l 



N I 



(43) 



(44) 
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Parameter 



LJ units 



Atomic Units 



tBB 
£AB 
OAA 
OBB 
GAB 

MassA 
Mass a 



1 

0.5 
1.5 
1 

0.88 
0.8 

1 

1 



3.8xl0 -4 
1.9xl0~ 4 
5.7xl0~ 4 
6.43 
5.65 
5.14 
3646 
3646 



Table I. Parameters used in our RPMD simulations on the 
Andersen-Kob Lennard- Jones glass forming system. 



The RPMD method has previously been used to study 
a diverse selection of multidimensional systems includ- 
ing proton transfer between organic molecules,^ diffu- 
sion in and inelastic neutron scattering from liquid para- 
hydrogenj 43 ! 44 ! diffusion of light atoms in liquid waterj^ 
and gas phase reactions such as that between methane 
and hydrogenPsl In all cases RPMD has been able to 
capture the dominant quantum mechanical effects in the 
dynamics and provide good agreement with the avail- 
able experimental or exact results. RPMD has also been 
applied to look at deep tunneling of Muonium and Hy- 
drogen atoms in ice^ and in this regime has been shown 
to be related to semi-classical Instanton theoryP^ 



V. SIMULATIONS DETAILS 

The quantum mode coupling theory requires as input 
the static structure factor and its Kubo transform. In the 
present study, we used a single component hard sphere 
(HS) model to generate this input within the frame work 
of the integral equation approach described above. Us- 
ing the PY closure, the system remains disordered even 
at very high volume fraction, thus providing a simple 
model to explore the quantum glass transition. The in- 
tegral equations (31l-p4|) were solved self-consistently. 



A simple trapezoidal integration scheme over the imag- 
inary time axis was employed, with P = 400 slices (we 
have checked convergence of the static input with respect 
to P). Here, P is analogous to the number of beads 
in the RPMD approach. For the HS system, it can be 
shown that the quantum mode coupling equations scale 
with the ratio of the de Broglie thermal wavelength to 
the particle size, A* = <J (3h 2 /ma 2 . Thus, to change the 
quantumness, one can either change h, or the mass, or 
the temperature. For the QMCT results shown below, 
we have varied the temperature to reflect a change in A* . 
We note that the temperature has no effect on the static 
structure factor in the classical case. 

We performed RPMD simulations on the Kob- 
Andersen glass forming systemj^liSI a binary LJ fluid, 
because the HS system investigated above by means of 
QMCT crystallizes on the timescale of the RPMD simu- 
lations. Each simulation consisted of 1000 particles, 800 
of type A and 200 of type B in a cubic box of length 



%Aoaa- The LJ parameters are given in Table |TjThe 
equations of motion were integrated using a time step of 
0.005 in Lennard- Jones (LJ) units using the normal mode 
integration scheme of Ref. 1501 The number of beads, P, 
used was given by the formula, 



P = 



11.2ft 



(45) 



This choice gives good convergence for all the regimes 
studied. Initial configurations were generated by anneal- 
ing from a temperature T* — 5.0 to the target temper- 
ature over a period of 2 x 10 6 time-steps. From these 
initial configurations we ran a further 2 x 10 5 steps of 
equilibration using a targeted Langevin equation normal 
mode thermostatting scheme This was followed by mi- 
crocanonical dynamics for 2 x 10 6 steps during which the 
results were collected. The quantum effect, A*, was var- 
ied by changing the parameter h. Five simulations were 
run for each temperature and value of h and the results 
averaged. 



VI. RESULTS 

Fig. [I] shows the results obtained from our QMCT 
treatment of hard spheres and RPMD simulations of the 
KA binary LJ fluid as the size of quantum fluctuations in 
the system are variedP Both of these systems have pre- 
viously been shown to exhibit all of the features of glassy 
behavior present in more complex fluids. In panel (b.) 
we show the liquid-glass dynamic phase diagram that is 
obtained from the QMCT calculation. The phase bound- 
ary is defined as the point where the solution of equations 
Eqs. (21 1, (22 1 and ( 29 1 leads to a finite value for the non- 
ergodic parameters, f q . At this point QMCT predicts 



that the system will never fully relax on any time-scale 
at the given packing fraction. For the RPMD calcula- 
tions, which are based on the evolution of semi-classical 
trajectories, we instead show the effect of quantum fluc- 
tuations on the diffusion coefficient of the particles at two 
different temperatures (T* =2.0 and 0.7) as the classical 
glass transition temperature of the system is approached 
(T* « 0.45) in panel (a.) of Fig. [T] Since the mean 
square displacement of the particles in the ring polymer 
trajectories show a caging regime (see the panel (c.) of 
Fig. [TJ, the diffusion constant was extracted from the 
long time slope of the mean-square displacement where 
the diffusive regime had been reached. The size of the 
quantum fluctuations were controlled by varying A* , the 
ratio of the de Broglie thermal wavelength to the particle 
size which controls the scale of quantum behavior. 

Comparing the RPMD results in panel (a.) and QMCT 
results in panel (b.) of Fig. [T| a remarkably consistent 
picture emerges from these two different approaches to 
quantum dynamics and glass forming systems. In the 
classical limit (A* — > 0) RPMD reduces to classical me- 
chanics and QMCT to classical MCT. As small quan- 
tum fluctuations are initially introduced, little difference 
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Figure 1. Panel (a.): The diffusion constant of particles of type A as a function of the quantumness, A*, obtained from 
the RPMD simulations for a quantum Kob- Anderson LJ binary mixture for two temperatures. Panel (b.): Dynamic phase 
diagram (volume fraction versus quantumness) calculated from the QMCT for a hard-sphere fluid. Panel (a): The mean square 
displacement of A particles as obtained from the RPMD simulations for the classical case (left frame, A* =0), the trapped 
regime (middle frame, A* = 1.125), and the regime governed by strong quantum fluctuations (right frame, A* = 1.1325). 



is observed in either the RPMD diffusion coefficient or 
QMCT liquid-glass line. However, as A* is increased be- 
yond 0.1, quantum effects are at first found to promote 
and then inhibit glass formation. In the case of RPMD, 
this is characterized by a decrease of nearly three orders 
of magnitude in the diffusion coefficient, and for QMCT, 
a 20 % fall in the packing fraction required for vitrifi- 
cation. When the thermal wavelength is increased fur- 
ther and becomes on the order of the particle size, the 
diffusion coefficient in the quantum system exceeds that 
observed in the classical limit. In addition the RPMD 
simulations at T* = 0.7 and 2.0 indicate that size of the 
re-entrance becomes much larger as the glass transition 
temperature is approached. Moreover, there is a hint of 
an interesting effect where, at high values of A*, the diffu- 
sion coefficient at lower temperature exceeds that at the 
higher temperature. We will return to this point later. 

Since both MCT and our new QMCT approach use 



the structure factor as input it is instructive to see if 
the dynamical reentrance is hinted at in this property. 
Fig. [2] shows the radial distribution function (RDF), 
which is the spatial Fourier transform of the structure 
factor, that has been calculated from the RPMD simula- 
tions of the KALJ fluid. For static equilibrium properties 
such as the RDF, RPMD gives numerically exact results 
since it reduces to the path integral molecular dynam- 
ics approach.^ The true (observable) quantum RDF is 
determined by the ring polymer bead correlations and 
is shown in the top panel for both the classical limit 
(A* = 0) and for a trapped regime (A* = 0.75). As quan- 
tum effects are introduced the RDF exhibits a broadening 
of the peaks due to the increasing uncertainty in the par- 
ticle positions which acts to smear out the pair structure. 
Throughout the entire range of A* studied the structure 
is observed simply to broaden systematically with A* and 
thus, there is no indication of the observed dynamical 
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Figure 3. Root-mean-square of the radius of gyration of A 
particles as a function of A* obtained from the RPMD sim- 
ulations for a quantum Kob- Anderson LJ binary mixture for 
two temperatures. The radius of gyration is defined as the 
average distance of the replicas from the polymer center. The 
results are plotted for temperatures T* = 0.7 (circles with 
dashed lines) and T* = 2.0 (triangles with dotted lines). 



Figure 2. The bead (upper panel) and centroid (lower panel) 
radial distribution functions of A particles for a classical 
(A* = 0, dashed) and trapped quantum (A* = 0.75, solid) 
regime. The bead distribution suggests less order in the 
trapped regime compared to a classical simulation while the 
centroid structure shows an increase in order. 



reentrance in the RDF. 

In the bottom panel we show the centroid RDF in 
which the centers of the imaginary time paths, rather 
than the bead positions, were used to compute the RDF. 
In the classical limit all beads collapse to a single point 
and hence both ways of calculating the RDF become 
identical. However as quantum fluctuations are increased 
the beads spread further from the center of the polymer 
and hence the centroid structure offers a different view 
into the structure of the quantum liquid. Upon examin- 
ing the centroid RDF in Fig |2]one sees the opposite trend 
upon increasing quantum fluctuations to that observed in 
the bead RDF, i.e. weak quantum fluctuations lead to a 
more structured centroid RDF which one would associate 
with more glassy dynamics. As quantum fluctuations fur- 
ther increase this trend reverses (data not shown) . Hence 
the centroid pair distribution function, which is not an 
experimental observable, appears to grossly mimic the 
dynamical correlations observed in both the QMCT and 
RPMD calculations. This is not entirely surprising, be- 
cause one expects that the centroid molecular dynamics 
(CMD) method, 1 ^ an approach similar to RPMD, will 
also capture the reentrance. Since CMD is an effective 
classical dynamics on the many-body centroid potential 



and since there are situations where the many-body cen- 
troid potential can be approximat ed by a sum of pair- wise 
potentials given by — fc^Tlog g(r)p^Jsuch static corre- 
lations in the centroid RDF must be evident if CMD 
is to reproduce the same phenomenology as predicted by 
QMCT and RPMD. This fact suggest that a strictly clas- 
sical MCT calculation that uses a static structure factor 
constructed from the centroid correlations might be a 
good proxy for the full QMCT calculation. It should be 
noted that the full QMCT only uses the observable struc- 
ture factor and thus one role played by the quantum ver- 
tex function is to effectively convert the bead correlations 
to centroid ones via the quantum fluctuation-dissipation 
theorem. The fact that the quantum vertex involves fre- 
quency convolutions while the classical version does not 
suggests, however, that there must be some distinction 
between a classical MCT calculation with centroid corre- 
lations and the full QMCT. 

So what is the origin of the reentrance? For this we 
turn to the RPMD trajectories to provide a physically 
insightful picture. Since this approximation maps the 
dynamics of a quantum mechanical particles onto that of 
a system of classical ring polymers, we can initially inter- 
pret the results in the language of the diffusion of classical 
polymers. In doing so we are careful to note that each 
bead on a given polymer only interacts with the bead on 
another ring polymer corresponding to the same imagi- 
nary time slice, a point which we will return to later in 
this section. In the non-interacting limit, the free ring 
polymer radius of gyration is directly proportional to the 
thermal deBroglie wavelength of the quantum particle. 
Hence, increasing A* allows the ring polymer represent- 
ing each quantum particle to spread out. The average 
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radius of gyration of each quantum particle in the inter- 
acting KALJ system is a static property which can be 
calculated exactly from RPMD simulations. In Fig.[3]we 
plot the average radius of gyration of each ring polymer 
relative to the value in the free limit. The dependence of 
this ratio on A* mimics the dependence of the diffusion 
coefficient on A* shown in Fig. [T] The decrease in this 
ratio when reentrance is observed suggests a correlation 
between the localization of the quantum particle and the 
increase in the glassiness of the system. As quantum fluc- 
tuations are increased from A* < 0.1, the effective diam- 
eter of the quantum particles differ little from a so that 
they can still fit into the thermally accessible space, their 
radius of gyration is still well approximated by B% ee , and 
little change in the dynamics is observed. However, once 
A* exceeds 0.1 there is not enough free space for the free 
ring polymers to further expand and crowding due to the 
surrounding solvent cage causes the radius of gyration to 
decrease from its free particle value. 

In the upper panel of Fig. [4j we show typical config- 
urations of a RPMD trajectory in the regime where the 
particle is localized in a cavity. The particle is confined 
by its surrounding neighbors, thus giving rise to an in- 
crease in its quantum kinetic energy. For diffusion to 
occur, particles must push past each other, causing fur- 
ther localization and incurring an even greater increase 
in their kinetic energies. This energy penalty to motion 
leads to slower dynamics. As A* is further increased, 
a tipping point is reached when the thermal wavelength 
becomes comparable to the particle size, A* « 1, At this 
point, the cost of localization becomes so large that the 
induced quantum kinetic energy enables the crossing of 
barriers between cavities, leading to a rise in the radius 
of gyration and facilitating diffusion. This can be seen 
in the representative snapshots of a RPMD trajectory 
shown in the lower panels of Fig. [4] in which the particle 
is delocalized across two cavities. Accordingly, the radius 
of gyration recovers with a corresponding increase in dif- 
fusion coefficient and diminishing of the caging regime. 
This can be likened to a "lakes to oceans" percolation 
transition, in which the caging regime reflects frustra- 
tion of the quantum particle in the classical potential, a 
frustration which is reduced when the kinetic energy of 
confinement essentially floods the barriers and allows the 
particle to traverse the region between adjacent potential 
energy minima. 

Reentrant effects in quantum systems have also been 
observed in the diffusion of electrons in a sea of classi- 
cal random blockers^ as well as in model systems In 
the former case the problem can be exactly mapped onto 
the diffusion of a classical ring polymer. However, in our 
case, while the expression "ring polymer" is used to de- 
scribe the isomorphism arising from the imaginary time 
path integral representation described in Eqs. (39l-(42l, 



it is not simply that of a system of true harmonic ring 
polymers. This is because each bead of a polymer only 
interacts with its corresponding bead at the same imag- 
inary time on the polymer representing another particle 




Figure 4. A series of snapshots taken from simulations at 
A* = 1.125 (upper panels) and A* = 1.3125 (lower panels) 
with T* — 0.7. For clarity the full imaginary time path (col- 
ored red) is only shown for one particle of type A with all 
others represented by their centroids. The centroids for the 
other particles of types A and B are colored green and blue, 
respectively. The upper panels depict configurations which 
reside in the trapped regime where the ring polymer is essen- 
tially localized in one cavity cage whereas in the tunneling 
regime (lower panels) it is frequently spread across two or 
more cavities in the liquid resulting in more facile motion. 



and not with any other beads on that particle. One might 
therefore expect that in systems with strong interactions 
it might be advantageous for the polymers to correlate 
their beads so as to minimize repulsion in exchange for 
a loss in entropy. To investigate this, we define vectors 
R„ = r\t^ — r„ , which represent the position of the bead 



,fc _ _(fe) 

at imaginary time k on ring polymer a relative to the 



position of the centroid (r° = (1/P) J2k=i r a )> an d we 
define the angle between vectors and R?) as, 



cos 9 



(fe) 

a,/3 



R 



r: 



K /3 



(46) 



(fe) 

This function, cos 9 a p, will have a value of —1 if the k-th. 
beads on polymers a and (3 are aligned perfectly away 
from each other and +1 if the beads are aligned towards 
each other. Since any correlation between the beads on 
two different particles is likely to be more pronounced at 
short distances where interactions are stronger we plot 
the correlation function C(r), 

c w - (jj E j> E cos - K P \)h ( 4? ) 

a>P k=l 

as a function of the distance between the centroids of two 
ring polymers. 

The function C{r) is shown in Fig. [5] for A* = 0.75, 
which corresponds to the trapped regime, and for A = 
1.3125, which corresponds to the strong quantum fluctua- 
tion regime. For r < a C(r) is negative in both cases. At 
these distances the potential between particles is strongly 
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Figure 5. The bead vector correlation (see Eq. (47l) for a 
trapped regime with A* = 0.75 (left panel) and regime where 
quantum fluctuations are pronounced with A* = 1.3125 (right 
panel). The solid lines represent the bead vector correlations 
between A particles and the dashed ones those between B 
particles. In both cases T* = 0.7. In the trapped regime the 
ring polymer beads show a large positive correlation around 
r — a which results in a large repulsion when the particles 
attempt to move past each other. In the other regime the 
beads align such that the correlation is largely negative which 
facilitates particle motion. 



repulsive and hence for polymers to approach this close 
their beads for the same imaginary time must avoid each 
other. However for r ss er, C(r) corresponding to the 
lower value of A*, the correlation becomes positive. At 
this distance the pair potential is attractive and hence the 
energy of the system is lowered if the polymer arranges its 
beads such that they are aligned on the same side of the 
respective ring-polymers. However at the higher value of 
A* . the entropic cost of such an ordering outweighs the 
energetic benefit, and hence C(r) is negative. This co- 
incides with the change between the dynamical regimes 
of quantum trapping and strong fluctuations because, for 
diffusion to occur, particles must move past each other. 
This regime corresponds to enhanced tunneling. In the 
case of low A* the beads of the polymer in the first co- 
ordination shell at r — a are largely aligned such that 
pushing them together induces a larger repulsion than if 
no such correlation existed. This increases the barrier to 
diffusion in this regime. 

One natural question that arises from this interpreta- 
tion of our results is what occurs if the RPMD calcu- 
lations are carried out at constant pressure rather than 
constant volume. The analogous QMCT calculations are 
constant volume calculations and indeed, as far as the 
hard-sphere control variables of volume fraction and A* 
are concerned, this question is irrelevant. The pressure 
varies as the volume fraction, which can then just be 
rescaled to yield results identical to those presented in 
the panel (b.) of Fig. [T] However from the standpoint of 
thermal variation, e.g. the variation of diffusion at fixed 
temperature while varying A* (see panel (a.) of Fig. [T]), 
this question needs to be addressed. A natural expecta- 



tion is that at constant pressure the reentrant effect will 
be mitigated or destroyed as the system can now adjust 
its volume as a natural response to the buildup of local 
pressure created by the "swelling" of the ring polymer. 
However, some aspects of this effect have been observed, 
for example, in analogous reentrant-like effects seen in 
Ref. 1451 where quantization of a single species in a con- 
stant pressure classical bath produces a reduction of the 
effective diffusion constant. More generally the values 
of A* for which the slowing of the liquid is observed are 
highly realizable in room temperature systems. The ther- 
mal wavelength of hydrogen at 300 K is 1.0 and hence 
the region of quantum slowing corresponds to diffusion 
in a medium with particles of radius 2 to 5 . Such a slow 
down is evident in experimental measurements of the dif- 
fusion of h ydrog en in non-glassy media such as water and 
palladium.HsHfiT 

We have carried out RPMD simulations of the binary 
glass-forming system at constant pressure, and indeed 
found at least a strong mitigation of the reentrant effect. 
Currently our statistics are not sufficient to make defini- 
tive statements about dynamical behavior in these sys- 
tems, and thus these results will be reported in a future 
publication. Regardless, it is clear that constant volume 
(confined) systems will exhibit a strong enhancement of 
the effects reported here. Further it should be mentioned 
that a similar reentrance is seen in lattice models of quan- 
tum glasses where the concept of swelling of imaginary 
time path s cannot be invoked to explain reentrant relax- 
ationPni 

On a final note, a subtle feature of the RPMD results 
of Fig. [I] (panel (a.)) should be mentioned. At very 
large values of A* the isothermal diffusion curves appear 
to cross. While the effect is quite small, this crossing 
would imply a reentrance of a different sort, namely a 
"melting by cooling" mechanism. This type of reentrance, 
distinct from that discussed for the bulk of this work, is 
similar to that discussed in Ref. §5J It should be noted, 
however, that the A* values here are large enough that 
particle statistics cannot be neglected in the simulation 
of a realistic quantum fluid, and the inclusion of such 
features may obviate this effect. 



VII. CONCLUDING REMARKS 

In this work we have presented a self-contained discus- 
sion of predictions for quantum glasses made by QMCT 
and RPMD. The predictions of these two distinct, albeit 
highly approximate, theories appear to be in harmony 
with each other. Both predict a strong reentrance in the 
relaxation of quantum supercooled liquids, namely that 
weak quantum fluctuations actually serve to push the sys- 
tem closer towards the glass transition. This seemingly 
paradoxical effect has also been noted in lattice models of 
quantum glasses and in models of quantum optimization. 
Indeed, one interesting aspect of our work is that it sug- 
gests that typical quantum annealing protocols should 
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generically have regions of parameter space where they 
are in fact less efficient than their classical counterparts. 

Future work will be directed towards the inclusion of 
bosonic statistics into the formulation of QMCT so that 
an investigation of the putative superglass may be car- 
ried out in a microscopic manner. In addition, it would 
be interesting to investigate more complex liquids such 
as confined supercooled water to see if quantum effects 
which may manifest at high temperatures lead to novel 
dynamical relaxation patterns. These topics will be re- 
served for the future. 
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